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Abstract 

We identify the global symmetries of SU(2) lattice gauge theory with N flavors of 
staggered fermion in the presence of a quark chemical potential \i, for fermions in both 
fundamental and adjoint representations, and anticipate likely patterns of symmetry 
j> \ breaking at both low and high densities. Results from numerical simulations of the 

model with N = 1 adjoint flavor on a 4 3 x 8 lattice are presented, using both hybrid 
Monte Carlo and Two-Step Multi-Boson algorithms. It is shown that the sign of the 
fermion determinant starts to fluctuate once the model enters a phase with non-zero 
baryon charge density. HMC simulations are not ergodic in this regime, but TSMB 
simulations retain ergodicity even in the dense phase, and in addition appear to show 
superior decorrelation. The HMC results for the equation of state and the pion mass 
show good quantitative agreement with the predictions of chiral perturbation theory, 
which should hold only for N > 2. The TSMB results incorporating the sign of 
the determinant support a delayed onset transition, consistent with the pattern of 
symmetry breaking expected for N = 1. 
PACS: ll.10.Kk, 11.30.Fs, 11.15.Ha, 21.65.+f 
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1 Introduction 



In QCD it is believed that as the baryon density rises the degrees of freedom providing 
the most suitable description change from being "hadronic" , ie. composite states such 
as protons and neutrons, to being "partonic", ie. quarks and gluons. This change may 
well be signalled by a phase transition as the appropriate thermodynamic variable, the 
baryon chemical potential //, is raised. Recent theoretical speculation |TJ, |2|, |3| suggests 
that the ground state of strongly-interacting quark matter at high density may be 
more exotic than initially thought, for instance existing in a superconducting and/or 
superffuid state due to the condensation of diquark pairs at the Fermi surface via a 
BCS instability. In such a state the SU(3) color gauge group is spontaneously broken 
by a dynamical Higgs mechanism; in the language of condensed-matter physics this is 
the Meissner effect. As well as being of intrinsic theoretical interest, the behaviour of 
strongly-interacting matter at extreme densities is of fundamental importance both to 
nuclear physics and in understanding compact astrophysical objects such as neutron 
stars. 

Unfortunately, it is difficult to apply the most reliable calculational tool for QCD, 
lattice gauge theory, directly to this problem. The reason is that once \x 7^ 0, the 
anti-hermitian property of the Euclidian Dirac operator D governing the motion of 
the quarks q and anti-quarks q in the presence of the color field is spoiled, with the 
result that the functional measure detM, with M = D([x) + m where m is the quark 
mass, is no longer positive definite. To apply the normal Monte Carlo method of 
importance sampling the path integral, the determinant must be split into a modulus 
and a phase; importance sampling is then done with respect to a measure |detM|, 
and arg(detM) is incorporated with the observable: 

((Oarg(detM))) 
{ ' ((arg(detM))) ' 1 ' ' 

where ((. . .)) denotes the expectation value with respect to the positive real measure. 
Now, the denominator of ( |1.1| ) is in effect a ratio of the partition functions of two 
different theories, one the true theory and the other an artificial one with a positive 
real measure: it should thus scale as exp(— AF), where the free energy difference 
between the theories AF is an extensive quantity. The number of states to be sampled 
before estimators of observables converge therefore in general rises exponentially with 
the system volume. This is the origin of the notorious 'sign problem'. 

We can gain physical insight into the sign problem by considering standard QCD 
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simulation algorithms such as the hybrid Monte Carlo (HMC) algorithm, which use a 
positive definite measure detM'M. The has the effect of introducing "conjugate 
quarks" q c transforming in the conjugate representation of the color group M. In 
general the presence of light gauge invariant bound qq c states carrying net baryon 
number in the spectrum results in unphysical behaviour for /i 7^ 0, eg. a premature 
"onset" transition between the vacuum and nuclear matter at /i D ~ 0{m 7T ) 1 the pion 
mass, which if chiral symmetry is spontaneously broken is much smaller than the 
constituent quark mass scale at which the transition is expected. There are, however, 
two strongly-interacting model theories where conjugate quarks can be tolerated. 
First consider the Nambu-Jona-Lasinio (NJL) model, long used as an effective theory 
for strong interactions; for sufficiently strong coupling it displays chiral symmetry 
breaking, signalled by a non-vanishing condensate (qq) and the development of a 
constituent quark mass greatly exceeding the current quark mass m, together with 
a triplet of light mesonic (ie. qq) pion states via Goldstone's theorem. Conjugate 
quarks in numerical simulations of the NJL model are harmless because the diagrams 
responsible for the tight binding and small mass of the pions are only accessible in qq 
channels ||. 

Next, consider non-abelian gauge theory in which the color group is SU(2) rather 
than the physical SU(3). Over the years "Two Color QCD" has been used to study the 
strong interaction in a variety of different contexts, the main motivation being that 
the computer effort required is appreciably less. Once quarks are introduced, however, 
important physical differences become apparent. Because the matter representations 
are all either real or pseudoreal, there is no gauge quantum number distinguishing q 
from q, resulting in enhanced global flavor symmetries in which both qq meson and 
qq baryon states appear in the same multiplet 0, and chiral (qq) and diquark (qq) 
condensates are related by global rotations [0, f| |§. The lightest baryon is therefore 
degenerate with the pion, and the onset chemical potential vanishes as \x Q oc y/m in 
the chiral limit. The same features ensure detM(/x) is real, and hence the theory 
simulable using standard algorithms || OH- 

Our motivations for studying Two Color QCD with adjoint quarks are twofold. 
Firstly, as we shall demonstrate, the pattern of symmetry breaking anticipated for 
gauge theories with quarks in real or pseudoreal representations of the gauge group 
differs between continuum || and staggered lattice fermions ||. There are good 
reasons, therefore, for considering the model with adjoint rather than fundamental 
quarks to be the most 'QCD-like' of the lattice models; in particular it has gauge 
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invariant spin-| states in its spectrum, for N — 1 quark flavors no baryonic Goldstone 
modes are expected, and it could potentially have a superconducting ground state at 
high density. Secondly, for an odd number of flavors the functional measure, though 
real, is not positive definite. The model thus has a potential sign problem of a simpler 
form than QCD in that only sectors having opposite sign, rather than a continuum of 
phases, need be considered. It may thus be feasible to make progress using standard 
means, or at least expose physical distinctions between the two sectors. 

Apart from the sign problem, which may almost be considered a problem of prin- 
ciple in the study of non-zero chemical potential, there are practical problems in using 
the HMC algorithm once \i ^ 0. For instance, the numerical effort required to invert 
the matrix M{fi) rises considerably as /i increases, due to the proliferation of complex 
eigenvalues with small modulus. This has encouraged us to study the performance of 
an alternative approach, the Two-Step Multi-Bosonic (TSMB) algorithm, in which 
M is not inverted at all, but the effect of detM incorporated by local Monte Carlo 
simulation over many auxiliary boson fields [|nj. Our conclusion is that the TSMB 
algorithm may be by far the more effective approach, in terms of the cost to produce 
decorrelated configurations, in the high density phase. Moreover for the class of prob- 
lem in which detM is real but not positive our results suggest the HMC algorithm 
in its simplest form is not ergodic in the same region of parameter space, since it 
fails to change the sign of the determinant. The TSMB algorithm does not share this 
problem. 

The remainder of the paper is organised as follows. In Section 2 we outline the 
global symmetries of the lattice model, and review the expected breaking pattern, for 
both adjoint and fundamental quarks. Operators for possible diquark condensates 
which form at high baryon density are discussed. We present a proof that detM is 
real and positive for fundamental lattice quarks, but only real for adjoint. Finally 
issues concerning the continuum limit, which is problematic, are discussed. In section 
3 we discuss both HMC and TSMB algorithms, and analyse their performance by 
considering the eigenvalue spectrum of M. It will be made clear that both sign 
problem and ergodicity are issues of practical importance once /i > /i , and that the 
TSMB approach is better suited to tackling them. In section 4 we present results 
of simulations using both algorithms. HMC simulations, confined to the sector of 
positive determinant, show clear evidence of an onset phase transition at fi Q ~ 
to a ground state with a non-zero density of baryon charge. The results moreover are 
in excellent agreement with analytic predictions obtained using chiral perturbation 
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theory ||, which might be expected to hold for more than one quark flavor. By way of 
contrast the TSMB simulations, which take the determinant sign into account, show 
evidence that this onset transition is thereby delayed. Our conclusions and future 
plans are briefly outlined in Section 5. 

2 Two Color QCD on the Lattice 

2.1 Formulation and Symmetries at Low Density 

In this section we review the formulation and symmetries of Two Color QCD with 
staggered lattice fermions in the presence of a chemical potential /x, and compare them 
to those of the corresponding continuum models. In this way we hope to motivate 
the study of the model with N = 1 flavor of adjoint quark as the most 'QCD-like' of 
the possibilities. The fermionic part of the lattice action is as follows: 

S = J2x p (x)D x AU,lAx p {y)+mr(x)6 x>yX p (y) = J2x P (x)M x>y [U,»)x P (y), (2.1) 

x,y x,y 

where the index p runs over N flavors of staggered quark, and D is given by 

+ \ Vo (x) (e»U (x)5 x>y _ t - e-»uZ(y)8 x>y+b ) . (2.2) 

The XiX are single spin component Grassmann objects, and the phases rj^x) are 
defined to be (-l) a o+-+^-i. 

In the case of fundamental quarks, the link matrices are complex 2x2 matrices 
acting on isodoublet x, x, and may be parametrised in terms of 3 real numbers aij as 
U = exp(iajTj), where 73 are the Pauli matrices. Note that T2UT2 = U*. For adjoint 
quarks, the same group elements may be represented by real 3x3 orthogonal matrices 
O acting on isotriplet %■> %■> given by 

On = hrinUr^). (2.3) 

In terms of the «j, O = exp(22«jf), where in this representation the generators 
{t l )jk — —i^ijk are hermitian, pure imaginary, and antisymmetric. For notational 
convenience we will continue to write the link variables as in either case. 

On integration over \ an d X the effective action exp(— S e ff) = det N M is obtained. 
Unlike three-color QCD with \x 7^ 0, exp(— S e ff) is real, since in the fundamental case 
detM = detr2Mr 2 = detM*, while in the adjoint case M is manifestly real. 
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In the chiral limit m — 0, the action has two manifest global symmetries: 



U(JV) e : X e 
U(JV) : Xo 



(2.4) 



the e/o subscripts denoting fields on even and odd sublattices respectively. However, 
it is straightforward to rearrange (|2.1|J2.2j) using the Grassmann nature of x-, X an d 
the fact that r]^(x ± jl) = Vni 30 ) to rewrite the action in this limit as 



1 



xeven,v 



X e {x) 
X e (x) 

where the fields X, X are given by 

X e = (Xe, -XeT 2 

for fundamental quarks and 

X e ^(x e ,Xe 



U u (x)X (x + i>)- 
Ut{x-v)X {x-v) 



X n 



x n 



Xo 

-T2X0 



Xo 
A o 



(2.5) 



(2.6) 



(2.7) 



in the adjoint case. In the limit \i — *■ the XJ(N) e ®XJ(N) symmetry thus enlarges to 
U(2iV): 

X a i-> VX X e i-> X e V ] V G U(2iV). (2.8) 

Note that the \J(2N) group emerges because exact symmetries are non-anomalous in 
lattice formulations; for a continuum model with Nf flavors the analogous symmetry 
enlargement is S\J(N f ) L ®S\J(N f ) R ®XJ(l) B -^SU(2N f ) §. 

For a non-abelian gauge theory we expect spontaneous chiral symmetry breaking 
to occur at low density, signalled by the appearance of a chiral condensate (xx) 0) 
which in the first instance we will consider to have the same form as the bare quark 
mass term in ( |2.1| ). To determine the pattern of symmetry breaking it is helpful to 
recast the condensate in terms of X, X: we find 



XX ^ 



x P 



t 



T2K + xi 



tr 



in the fundamental case, and 



XX 



x P 



XT - x. 



tr 



r 2 X 



X, 



(2.9) 



(2.10) 
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for the adjoint. Here 1 denotes the N x N unit matrix. The residual symmetry left 
unbroken by the condensate is that which leaves invariant respectively the symmet- 
ric/antisymmetric 2N x 2N form. For \x = we thus find 

fundamental : U(2iV) -> 0(2N) adjoint : U(2iV) -> Sp(2iV). (2.11) 

This is remarkable in that it is the opposite of the breakdown in the continuum |J: 

fundamental : SU(2N f ) -> Sp(2N f ) adjoint : SU(2iV» -> 0(2iV/). (2.12) 

In effect the roles of fundamental and adjoint representations are reversed for stag- 



gered lattice fermions, a fact which has been noted several times over the years [jl2 

Next we identify the massless modes arising in the m, \x — > limit as a result of 
Goldstone's theorem. For the fundamental case, the number of broken generators 
predicted by (|2.11|) is N(2N + 1). The case of N — 1 fundamental quark has been 
analysed in [|J by considering infinitesimal rotations of the condensate fl2.9|) by V$ = 
1 + iS\, with A one of the U(2) generators {1,7;}, and identifying the Goldstone 
mode with the coefficient of 0(5). The three Goldstones thus found are the familiar 
mesonic 0~ pion x £ X (where e(x) = (— l) x o+ x i+ x 2+ x z^ an d two scalar + diquark 
states x ir T 2 x, X T 2X tr - The ± superscripts here denote the symmetry of the state 
under the following lattice 'parity' symmetry: 

x = (x ,xi,x 2 ,x 3 ) i-> x' = (x , 1 -Xi, 1 -x 2 , 1 -x 3 ) 
X (x) » (-l^xM 5 X(x) ~ (-lf^x(x'). (2.13) 

For m ytz o, the three states remain degenerate ||, gaining masses oc y/m in 
accordance with standard PCAC arguments. In the case of adjoint quarks, the pattern 
( P- 11| ) predicts N(2N — 1) Goldstones in general, and for N — 1 the only Goldstone 
mode is the 0~ mesonic pion. The corresponding analysis for the continuum models 
has been performed for Nf > 2 with some care in ; the Goldstone counts are found 
to be N f (2N f - 1) - 1 (fundamental) and N f (2N f + 1) - 1 (adjoint). Modulo the 
mode destroyed by the U(l) axial anomaly, the reversal of fundamental and adjoint 
cases is clear. Another important distinction is that in the continuum the Goldstone 
spectrum always contains diquark states; in this respect the N — 1 lattice adjoint 
model is special. 

Now consider the effect of increasing \x from zero. The chemical potential has the 
effect of promoting a ground state containing baryonic matter, signalled by a non-zero 
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value for the baryon number density 

n = \ (x(x)7] (x)[e^U (x)x(x + 0) + e -»U}(x - 0) X (x - 6)]) . (2.14) 

At zero temperature n thus serves as an order parameter for an onset phase transition 
occuring at some \x a separating the vacuum from a state containing matter. A naive 
energetic argument would suggest that the onset transition should occur for a value 
of /j, equal to the mass per baryon charge of the lightest particle carrying non-zero 
baryon number. For the models discussed in the previous paragraph in which some of 
the Goldstone modes are diquark states, those states will be the lightest baryons in 
the spectrum. Generically then, we expect /i D ~ m n /2 for most variants of Two Color 
QCD, which should be contrasted with the much larger value m nuc i eon /3 expected in 
physical QCD. The exception is the lattice adjoint model with N = 1. The fact that 
the lightest baryons are Goldstone modes and thus amenable to analysis by using 
Chiral Perturbation Theory (xPT) on an effective sigma model has been exploited in 



13| to calculate both Goldstone spectrum and equation of state of the continuum 



models as functions of //; the principal result is the prediction /x G = m n /2. 



2.2 Diquark Condensation at High Density 

For a sufficiently high density of baryon charge, regardless of the nature of the bound 
states in the low energy spectrum, the fermionic nature of the quarks should ensure 
that the dominant degrees of freedom are governed by a Fermi-Dirac distribution. 
At zero temperature, all states will be occupied up to the Fermi energy Ep, which 
coincides with [i in the limit where inter-quark interactions can be neglected (for an 
asymptotically free theory we expect this approximation to improve as /x rises). The 
question now arises as to whether this simple description is unstable with respect 
to condensation of diquark pairs situated at antipodal points on the Fermi surface, 
resulting in an energy gap between the ground state and the lowest spin-| excitation. 
For large \i this instability is generic provided the quark-quark interaction is attrac- 
tive; that this is so for non-abelian gauge theories has been argued as arising from 
either gluon exchange [I], or instanton effects Physically, diquark condensation 
implies that fermion number and/or baryon charge is no longer conserved, and the 
ground state is a superfluid. For physical QCD, diquark pairs cannot be color singlet, 
so the condensation results in the phenomenon of color superconductivity, rendering 
some or all of the gluons massive via a dynamical Higgs mechanism || [|. 
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We now consider possible diquark condensates that might form in Two Color 
QCD. There are many possible diquark states that can be written down, and in the 
absence of a detailed dynamical calculation we have to proceed by making some ad hoc 
assumptions. We might imagine that the wavefunction of the qq condensate should 
ideally be gauge invariant, spacetime scalar, and in the case of the lattice model, as 
local as possible in the x fields, since non-local wavefunctions require the insertion 
of link variables to maintain gauge invariance, whose fluctuations will weaken the 
condensation. The most important consideration (and the only one which is inviolate 
0, [TJ]]) is that the condensate respects the Pauli Exclusion Principle, implying that 



the operator be antisymmetric with respect to exchange of quantum numbers between 
the quarks. 

For Two Color Lattice QCD with fundamental quarks, a diquark operator which 
satisfies all of the above requirements is 

QQ2 = 2 [x U \x)t 2 x{x) + x(x)r 2 x tr (x)] . (2.15) 

In fact ||, qq2 is related to the chiral condensate XX by the U(2iV) symmetry (|2.8| ). 
and as \x increases the chiral condensate ( |2.9| ) is in effect rotated into the diquark one. 
There is, however, a physical distinction between high and low density phases; with 
N set to 1 and m, \i > the U(2) symmetry is reduced to U(l)s which is associated 
with conservation of baryon number. Condensation of qq-i spontaneously breaks this 
residual symmetry resulting in an exactly massless Goldstone mode. Because of the 
change in the number of massless modes, there is a true phase separation between 
either the vacuum or a low density normal phase, and a high density superfluid phase. 
Preliminary Monte Carlo simulations have revealed a transition as fi is increased 
f|, and the associated condensation of qq^ has been confirmed fl5|]. A superfluid 
condensate at high density is also found in both fundamental and adjoint continuum 
models ||; here the sigma model approach permits a calculation of the spectrum and 
equation of state as a function of /i/m, and predicts that (qq) and n become non-zero 
at the same critical \x Q = m n /2, ie. there is no normal phase. Another important result 
is that for \i ^> m the pseudo-Goldstone states, which would have been massless in 
the SU(2Aj) symmetric limit, have mass 2/i P, [13 . 



Next we discuss the possibilities for Two Color Lattice QCD with adjoint quarks 
|l6fl . For N > 2, a diquark operator satisfying all of our ideal requirements can always 
be written down: for N = 2 it reads 

i 



X ptr (x)e^ X q (x) + x P (x)e Pq X qtr (x) , (2-16) 
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where p, q — 1,2 are explicit flavor indices. In fact, just as in the fundamental case, 
qqs is related to the adjoint chiral condensate (|2.10|) by a global symmetry, in this 
case the U(4) rotation given by 

V=l( i P p %) with J>=(J "J). (2.17) 

Condensation of qq 3 breaks U(l)# but leaves unbroken an SU(2) of isospin. Since 
there are diquark states among the 6 Goldstones expected from the breaking U(4) — > 
Sp(4), we expect the usual scenario to apply, with a transition to a superfluid phase 
at /i Q = m n /2. 

For N = 1 the Exclusion Principle prevents us from writing a diquark operator 
satisfying all the requirements, since X tr { x )x{ x ) = 0- We have considered two possi- 
bilities. Firstly, a non-local operator which is gauge invariant and scalar under Q2.13D 
is 

Ms = ^E^W(-l)^ \x tr (x)U,(x) X (x + fi)-x(x)U,(x)x tr (x + fi)} . (2.18) 

The (— 1) Xm factor ensures that qq' 3 is antisymmetric with respect to spatial exchange 
of x fields. In terms of the X, X fields, 

qq' 3 = lj2v^)(-l)^X e (x)U,(x)(^ _1 ) X G (x + £). (2.19) 



In the m, fi —>■ limit the global symmetries left unbroken by (qq' s ) ^ are a U(l) £ 
symmetry, 

X^e~ ia X ; X^Xe ia , (2.20) 
which is broken when m/0, and an 0(2) symmetry rotating x into x : 

X^QX ; X^XQ tr ■ Q=( COS ^ sin ^ . (2.21) 

\-sinf cos v J 

The pattern of symmetry breaking is thus U(2)— ►U(l)(g>U(l) e , leaving two unbroken 
generators and hence two Goldstone modes, which turn out to be parity even meson 
and diquark states, the diquark remaining massless once m, \l > 0. Since for m/0 
there are no exact Goldstones in the low density phase where chiral symmetry is 
broken, we once again predict a phase separation between a low density phase and a 
high density superfluid. Recalling, however, that the only Goldstone at low density 
is a qq meson, we are unable to apply the effective theory arguments of || and 
in this case do not expect \i Q = m^/2. 
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Secondly, consider an operator which is local but not gauge-invariant: 

1 



QQsc = 2 



X tr (x)t l X (x) + x(x)fx tr (x) , (2.22) 



whose consistency with the Exclusion Principle is now due to the antisymmetry of 
the generators t l . Under a gauge transformation, qq l sc transforms in the adjoint rep- 
resentation of SU(2), 

Qqlc ^ O tJ qql c , (2.23) 

which follows from the property U~ 1 \ t U = OyX* , true for arbitrary representations 
generated by A*. Therefore qq l sc acts like an adjoint Higgs field, and its condensation 
breaks the SU(2) color group to U(l), as in the Georgi-Glashow model of electroweak 
physics pTTfl . As the subscript implies, this is therefore a superconducting solution; 
ironically, in this case the superconducting phase is characterised (at least in per- 
turbation theory) by a massless photon. Once again, because of the change in the 
number of massless particles between the low density confined phase and the high den- 
sity superconducting phase, a true phase separation is expectedF], which is consistent 
with general properties of gauge theories with Higgs fields in the adjoint representa- 



tion |]19 |, and in possible contrast with the 'color-flavor locked' state anticipated in 
three-flavor QCD [|J, where continuity between high and low density phases has been 
postulated p0|. 



2.3 The Sign of the Determinant 

As discussed in subsection |2.1| , it is straightforward to show that the path integral 
measure, proportional to detM , is real. To determine whether it is positive definite 
requires a more subtle argument. Let us first make some general observations: 

• Lemma 1: let M be any diagonalisable operator and K be the complex conju- 
gation operator. If there exists a unitary operator T such that [KT, M] = 0, 
then detM is real. 

• Proof: let ip be an eigenvector of M with eigenvalue A; then ip = KTip is an 
eigenvector with eigenvalue A*: 

Mip = Xip > KTMip = KTXi/j > MKTijj = \*KTip M?p = A*^; 



1 Note that in the 2+1 dimensional SU(2) adjoint Higgs model, there appears to be no phase 
separation, since in principle the photon can acquire a mass via non-perturbative effects; it appears 
extremely light in the "Higgs" phase, however [jL8[. 
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hence detM is real, since the product over eigenvalues can be organised in 
complex conjugate pairs. 

Now this does not imply that detM is positive, since for real A, if> might be propor- 
tional to if>, and hence we cannot exclude the possible existence of non-degenerate 
real eigenvalues which may include an odd number of negative ones. To prove that 
if> and ip are linearly independent, we need KT to have another property. 

• Lemma 2: if (KT) 2 = —1, if) and if) are linearly independent. 

• Proof: 

(^) = (ifj\KTip) = (Tip\TKTip) = ((KT) 2 ip\KTip) = -(ip\i>) = 0, 

where we have used (<p\if>) = [Kip\K(f>). 

Hence if [KT,M] = and (KT) 2 = —1, all real eigenvalues are doubly degenerate, 
and detM > 0. 

The square of the anti-unitary operator KT determines the Dyson index (3 of M, 
the case (KT) 2 = 1 implying (3 = 1 and (KT) 2 = —1 implying /3 — 4 JJ. Physical 
QCD with three colors and fundamental quarks has no operator equivalent to KT, 
corresponding to /3 — 2. Now let us identify the operator T for both continuum and 
for staggered lattice fermions. In the continuum, 

M=(d u + igXAl)^ + /i 7o + m, (2.24) 

with A l u the gauge potential and X 1 the generator appropriate to the quark represen- 
tation. We find 

fundamental: T = C^ 5 ®r 2 , (KT) 2 = 1 (3 = 1; 

adjoint: T = C 75 <g> 1, (KT) 2 = -1 (3 = 4. (2.25) 

Here, C is the charge conjugation matrix defined by its operation on Euclidean her- 
mitian 7-matrices: C / ~f fJi C~ 1 = — 7^- For staggered lattice fermions M is given by 
QS> and 

fundamental: T = r 2 , (KT) 2 = -1 =>■ (3 = 4; 

adjoint: T = 1, (KT) 2 = 1 (3 = 1. (2.26) 

Once again, we see the roles of fundamental and adjoint representations reversed for 
staggered lattice fermions. We thus have a proof that the functional integral measure 
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is positive definite for continuum adjoint quarks and fundamental lattice quarks. 
There is no such proof for continuum fundamental quarks and lattice adjoint quarks, 
and as we shall demonstrate in section |3.1| , there are indeed isolated real eigenvalues 
and hence a sign problem for the adjoint lattice model at large chemical potential fi. 

2.4 The Continuum Limit 

In this paper our viewpoint will be to study the lattice model in its own right, that is 
as a strongly coupled theory with the potential to show superfluid or superconducting 
properties at high density, and not to attempt either continuum or chiral limits. We 
have therefore focussed on the global symmetries appropriate to staggered fermions. 
This should not obscure the fact, however, that there are interesting issues related to 
the continuum limit whose resolution is still not clear. The model with iV flavors of 
staggered fermion should correspond to a continuum theory with Nf = AN physical 
quark flavors. For the model with fundamental quarks, therefore, we expect the 
U(2AT) global symmetry to enlarge to SU(8iV), which is spontaneously broken by 
a chiral condensate to Sp(8AT): for the adjoint model the corresponding pattern is 
SU(8iV) — > 0(8 AT). Does the role reversal of fundamental and adjoint representations 
cease at some point? Since the behaviour of both is qualitatively similar in the 
continuum formulation, this is plausible [f| . On the other hand, we have argued that 
the lattice model with N = 1 staggered flavor exhibits a distinct behaviour, with no 
premature onset transition, and the possibility of a superconducting condensate. If 
this were the case, would the superconducting phase survive the continuum limit? 
Above we also argued that the sign of the fermion determinant detM is not positive 
in all cases, implying that it is important whether A" is even or odd, and once again 
continuum and lattice models have distinct behaviour. 

Let us outline a possible route to resolving these issues. First consider the effect 
of extending the continuum operator T of eqn. (|2.25|) by including the effects of a 
four-index flavor structure. One can then consider an operator 

T' = T ® T (2.27) 

where T is a 4 x 4 matrix acting on flavor; the case T = C75 has the property that 
(KT') 2 = —(KT) 2 , which would render the fermion determinant for 4 fundamental 
continuum flavors positive (as must clearly be the case). Next consider the form of 
the lattice T operator ( |2.26| ) in a basis of fields q aa , q aa carrying both a spinor index a 
and a flavor index a (implicit in the staggered fermion approach), each taking values 
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from 1 to 4 [21]. The appropriate structure is 



T q = C 75 ® C l5 



(2.28) 



where the first C75 acts on the spinor index, the second on flavor, and the final T2 (1) 
acts on color in the case of fundamental (adjoint) quarks. At non-zero lattice spacing, 



due to a term which is formally 0(a) when M is expressed in the q, q basis [21], this 
is the only exact symmetry of the form (|2.25|J2.26| ); as we saw, (KT q ) 2 = 1 for adjoint 



quarks, implying a possible sign problem. However, in the continuum limit we expect 
flavor symmetry to be restored, implying that T' q with the second C75 substituted by 
arbitrary V will also define symmetries of M. The existence of merely one T' q with 
(KT' q ) 2 = —1 (eg. T = 1) suffices to prove detM positive due to degeneracy of real 
eigenvalues; we expect as a consequence an enlargement of the Goldstone manifold to 
recover the continuum symmetry prediction. Away from the continuum limit, these 
extra baryonic Goldstone modes remain massive due to lattice artifacts. 



We can also examine the diquark condensates postulated in subsection |2.2| in terms 
of q, q | 14}j . In this representation the local diquark operators (|2.15| , [2.16| , |2.23|) take 
the form 

qq=\ \q tr (C l5 ® C75 ® T)q + q(Cj 5 ® C l5 ® T)q tr ] , (2.29) 



■ [q tr (C l5 ® C75 ® T)q + q(C lb ® C l5 <g> T)q tr \ , 

where the first operator in the tensor product in each term acts on spinor indices, 
the second on flavor, the third on gauge indices (or explicit flavor indices in the case 
of fl2.16| )), and all three are antisymmetric matrices. In the same basis the parity 
transformation ( 2.13|) reads 



q(x) ( 7o ® 75 ® l)q(x') ; q(x) i-> g(x')(7o ® 75 ® 1). 



(2.30) 



We thus confirm that the first C75 is consistent with the local condensates being 
spacetime scalars, the second C75 implying that the condensates transform in the an- 
tisymmetric tensor representation of the continuum flavor SU(4) symmetry expected 
away from the chiral limit, which has dimension 6. The situation is different for the 
non-local condensate ( [2.18 ), which reads 



Ms 



£ q tr (C ® C% ® l)q - q(C ® C 7 ! ® 1)? 



,-7tr 



(2.31) 



Invariance of qq' 3 under parity depends on the non-trivial action of (|2.30|) on flavor, 
and in fact the spacetime structure of the operator implies that this condensate is 
pseudoscalar. The symmetry of qq' 3 may be checked noting that C7* is a symmetric 
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matrix, in turn implying that this condensate transforms as a symmetric 10 of flavor 
SU(4). Without further detailed dynamical input, it is difficult to proceed; it may 
well prove that the route to the continuum is complicated, with several distinct phases 
found as the parameters (3, m and \i are tuned. 

Finally, we note that because adjoint sources screen color four times more effec- 
tively than fundamental ones, asymptotic freedom is lost in Two Color QCD with Nf 
massless adjoint quarks for Nf = AN > 11/4, ie. even for N = 1. Therefore in the 
chiral limit even the location of the continuum limit is a priori unknown; possible be- 
haviours include a continuum limit described by a non-perturbative renormalisation 
group fixed point at some finite value of the gauge coupling /3 C , or no interacting con- 
tinuum limit existing at all. In section |] we will demonstrate that these considerations 
are irrelevant at the values of j3 and m considered in our work. 

2.5 Summary 

In this section we have reviewed the formulation and symmetries of Two Color lat- 
tice QCD at length, and compared and contrasted it with corresponding continuum 
models. We conclude that the version of Two Color QCD which is most 'QCD-like' 
with respect to non-zero chemical potential fi and hence worthy of study, is that 
with N = 1 flavor of staggered lattice fermion in the adjoint representation, for the 
following reasons: 

• There are no Goldstone diquarks in the spectrum at zero density, and hence no 
reason to expect a premature onset transition at fi Q ~ m n /2. 

• The spectrum in the confined phase should contain gauge invariant fermionic 
bound states, either of a quark and a gluon, or of an odd number of quarks. 
Hence there is the possibility of a nuclear liquid phase, and the formation of 
a Fermi surface, before the restoration of chiral symmetry expected as fi is 
increased. 

• The path integral measure is real but not positive definite for any odd N; 
that of QCD is complex. It may prove possible to expose the importance of 
configurations with non-positive definite detM in the path integral with /x ^ 0. 

• There is a possibility of a gauge-variant diquark condensation at high densities, 
leading to color superconductivity. No other lattice model capable of being 



studied with /i ^ appears to share this feature M, 14, 15 . 
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It is conceivable that some or all of these factors are intimately linked. 



3 Simulation Algorithms 

We have studied Two Color lattice QCD with adjoint staggered fermions using two 



different simulation algorithms, the hybrid Monte Carlo (HMC) algorithm [ 22] , and a 
Two-Step Multi-Bosonic (TSMB) algorithm [23]. Here we outline the two methods, 
and point out some important features of each. 



3.1 The Hybrid Monte Carlo Algorithm 

The HMC algorithm starts from an expression for the action S in terms of bosonic 
auxiliary pseudofermion variables $: 

1 N 8 
S ^EE^W^O^.^fe) - f E TrtVO^ (3-1) 

^ x,y p=l ^ x,fi<v 

where the fermion matrix M is defined in fl2.lU2.2f ), (3 is the gauge coupling constant, 
and the trace over the plaquette is taken in the fundamental representation. In 
contrast to QCD, the $ may be taken real, since M is real. Gaussian integration 
over $, which is convergent since the eigenvalues of (M tr M)~ l are real and positive, 
yields a factor proportional to 

exp(-S eff ) = (VdetM tr M) N = (detM^. (3.2) 

This coincides with the correct functional measure det^M if detM is positive. In 
this case HMC correctly simulates N — 1 flavor of staggered lattice fermion. Note 
that for n ^ 0, the usual trick of evaluating (M tr M)~ l on just even lattice sites, thus 
reducing the effective number of fermion degrees of freedom by a further factor of 
two, is not available since the chemical potential spoils the anti-hermitian form of D, 
the off-diagonal part of M. 

It is instructive to analyse the performance of the algorithm by considering the 
eigenvalue spectrum of M. Because D only connects even sites with odd and vice- 
versa, if ij)(x) is an eigenvector with eigenvalue A = m + k, then so is e(x)ip(x), with 
eigenvalue A = m—n. Because D is a real matrix, if ip is an eigenvector with eigenvalue 
A then so is ip* with eigenvalue A*. Now for \x = 0, the matrix D is anti-hermitian, 
implying its eigenvalues are pure imaginary; we deduce that in this case the spectrum 
of M lies along the line A = m + iz, with z real. For /i ^ 0, the symmetries between A, 
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A* and 2m — A persist, but now the spectrum swells out to occupy a roughly elliptical 
region of the complex plane, whose horizontal dimension grows with increasing \i. 
The spectrum from a representative configuration is shown in Fig. [I]. A feature of 
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Figure 1: Eigenvalue spectrum from a 4 3 x 8 lattice with (3 = 2.0, m = 0.1 and 
/i = 0.35. Note the mismatch in scale between real and imaginary axes. 



interest is that even once the spectrum has swelled out, an appreciable fraction of the 
eigenvalues remain on the line A = m + iz; this contrasts with behaviour observed in 
three-color QCD, where all eigenvalues leave the line once /z 7^ but is similar 
to the spectrum of the random matrix Dirac operator with /i / and Dyson index 
(3=1 (corresponding to a matrix D with off-diagonal elements chosen from a chiral 
Gaussian Orthogonal Ensemble) P5| . Secondly, for this particular value of /i = 0.35 
the spectrum has broadened sufficiently for the extremal eigenvalues to have negative 
real parts. 

The HMC algorithm works by evolving the {U} fields in fictitious time r for 
a period called a trajectory, whereupon the resulting configuration is accepted or 
rejected via a Metropolis step; the acceptance probability is related to how well 
the pseudo-Hamiltonian flow conserves energy. At the start of each trajectory the 
pseudofermion and conjugate momentum fields are refreshed from a Gaussian heat- 
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Figure 2: Spectral flow under HMC evolution. 

bath, which ensures ergodicity in most cases. The trajectory length may be held 
constant or picked at random, but is usually chosen to be 0(1). We find, however, that 
once the eigenvalue occupation region has swelled to include the origin, the resulting 
presence of very small eigenvalues causes numerical instabilities, and necessitates 
smaller trajectory lengths — the data of Fig. [1| were generated using an average 
trajectory length of 0.18. An HMC update step acts on the configuration as a whole, 
and high acceptance is maintained by making a sequence of small changes: HMC is 
thus a small step-size algorithm. In Fig. |2| we show the spectral flow under HMC 
evolution in a small region close to the origin, starting from the configuration of 
Fig. [l] and evolving for five very short trajectories of length 0.003. For clarity only 
the eigenvalues from the first, third and fifth steps are shown. Note that there are 
a few eigenvalues which are real: a pair close to 0.1±0.08 on the first configuration, 
which increases to eight by the fifth. Close inspection reveals that between the first 
and third configurations a pair of eigenvalues jumps from the line A = m + iz at 
approximately 0.1 ± 0.03i to the real axis at ~ 0.1 ± 0.02 — such processes change 
the number of real eigenvalues by ±2, and the number on the line A = m + iz by 
=1=2. Between the third and fifth configurations, two pairs of conjugate eigenvalues 



18 



coalesce and then move independently along the real axis at ~ 0.1 ± 0.04 — this 
process changes the number of real eigenvalues by ±4. Finally between the third 
and fifth configurations we can also see pairs of eigenvalues coalesce and then move 
independently along m+iz at ~ 0.1±0.04i, changing the number of eigenvalues on this 
line by ±4. All other eigenvalues in the region, at generic points in the complex plane, 
merely exhibit small fluctuations under HMC evolution. The symmetries between A, 
A* and 2m — A are maintained at all times. 

Fig. H demonstrates unambiguously the existence of non-degenerate real eigenval- 
ues of A, and the fact that there is always an even number of them. There is no 
obstruction in principle to there being an odd number of negative real eigenvalues, 
resulting in detM < 0. However, since the processes discussed in the previous para- 
graph always create real eigenvalues in pairs at the same point, the only route to 
detM < can be if an isolated eigenvalue moves along the real axis and through 
the origin, forcing detM to evolve smoothly through zero. In the neighbourhood of 
such a point the effective action S e jf diverges, implying a strong repulsive force in 
the Hamiltonian flow, and hence a large kinetic barrier to changing the determinant's 
sign. This barrier is a feature of any model in which detM is real; in QCD, which has 
detM complex, it should be straightforward to change the sign by moving through a 
sequence of configurations in which the phase of the determinant gradually increases 
from to 7T. 

In Fig. |3| we show evolution of In detM for 40 representative configurations, each 
separated by 50 trajectories of average length 0.05. The simulation parameters have 
been chosen such that the left-hand edge of the spectrum extends well beyond the 
imaginary axis, so that we expect a non-vanishing probability of obtaining negative 
real eigenvalues. We see on inspection of the first, solid curve that HMC evolution 
does not appear to change the sign of detM. This could in principle be because the 
volume of configuration space with an odd number of real negative eigenvalues is 
very small; we can eliminate this possibility by considering a comparable sequence 
of quenched updates (dotted curve). In this case the sign of detM is seen to change 
frequently, confirming that such configurations exist and are relatively easy to find. 
Moreover, if one of the configurations with detM < is then reequilibrated with 
HMC and allowed to evolve, we see that once again the sign remains stable and 
negative (dot-dashed curve); this behaviour is observed to persist for both signs for 
approaching 20000 trajectories. We reach two important conclusions: 
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Figure 3: Evolution of |lndetM| x sign(detM) using HMC, quenched and TSMB 
updates on a 4 3 x 8 lattice at (3 = 2.3, m = 0.1 and fi = 0.6. 

• There are regions of parameter space with /i / for which detM can take 
negative values, ie. there is a sign problem. 

• The kinetic barrier at the origin prevents the HMC algorithm from changing 
the sign of detM, and therefore from exploring the whole of the system's config- 
uration space: in other words, the HMC algorithm is not ergodic in this region 
of parameter space. 



The sign problem is well-known [26] and can be addressed, at least in principle, by 
including the sign of detM with the observable as in this may be expected to be 
effective provided the average sign is significantly different from zero. The problem 
with ergodicity is less well-known — it can be anticipated in any model where M 
is real but its eigenvalues A complex; another example which would be interesting 



to study is the lattice Gross-Neveu model with discrete Z 2 chiral symmetry [27| . 
Note that it remains a problem of principle for even N, since we can still classify a 
configuration by the number of negative real eigenvalues of M. It is an interesting 
open question whether or not the absence of transitions between odd and even sectors 
is of practical importance. 
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Despite its problems, we have made an extensive study of Two Color QCD using 
the HMC algorithm, and the results will be surveyed in section |j. We next turn to 
an algorithm which has the potential to overcome the ergodicity problem. 



3.2 The Two-Step Multi-Bosonic Algorithm 



For the Two-Step Multi-Bosonic algorithm |23| we also need a hermitian fermion 
matrix. First one might consider 



QGu) = eM(ji) = Qi-fiY, 



(3.3) 



but for non-zero chemical potential this is still not hermitian. Hermiticity may be 
achieved by doubling the matrix size: 



Q{-fi) 
Q((j) 



Q(//)t . 



(3.4) 



Following from 



/ Q(-H)Q(H) 

V o 

/ M(//)tM(/i) 
I 



. 

Q(^)Q(-fi) J 


M{-n)^M{-n) 



and noting that 



(3.5) 



(3.6) 



det M(ji) = det M(/x)* = detM(-//), 

(where the first equality holds because M = D + m is real, the second because 
D tr (fi) = —D(—fi), and the spectrum of D is symmetric about zero), we deduce 



det Q(/i) 2 = {detM^u)} 4 , |detM(/i)| = {detQGu) 2 } : 



(3.7) 



In multi-bosonic representations of the fermion determinant [TT| polynomials are 
used to approximate the necessary inverse powers of x = Q(fi) 2 over some prescribed 
range x G [e, A]: 

1 



1 



£4 



P X) 



det M(/x)| 



detP(Q(//) 2 ) ' 



(3.8) 



This shows that in the present case the same polynomial approximations can be used 
as in recent numerical simulations of supersymmetric Yang Mills theory ^8 . 
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The usual way to represent the determinant of the polynomial in ( ft.8|) is by 
functional integration over complex pseudofermion boson fields. Since the fermion 
matrix is real, it is also possible to use a real multi-bosonic representation, just as in 
the HMC case. For this we need a different polynomial approximation: 

— ~P(x), |detM(/i)|~ r . (3.9) 



12 



{det P(Q(/x) 2 )}' 



Since the polynomial P(x) is supposed to have complex conjugate pairs of roots, one 
can decompose it, with an overall factor tq, as 

HQ^Y) = r T[[Q(j*) 2 + r 3 ] = r ]J[(Q^) + N f + u]} , (3.10) 

j 3 

with fij, Uj real. In order to achieve this form one has to choose the signs of the 
square roots of two complex conjugate roots appropriately: rj = {ijij + Vj) 2 , r J+ i = 
r* = (— ifij + Vj) 2 . The multi-bosonic representation with real pseudofermion fields is 
then 



1 

oc 



, / [*&] exp f - £ H tr [(Q(») + Mi) 2 + tfyH] ■ (3-11) 
{det P(Q(/i) 2 )} 2 J V iv* J 

The polynomial orders for a sufficiently good approximation in ( |3.9|) are typically 
somewhat higher than in ( |3.8|) but the use of real fields has the advantage of taking 
half the storage and roughly half the arithmetic. 

In the TSMB algorithm the polynomials in eqs. ( |3.8|) and (|3.9| ) are obtained by 
a product of lower order polynomials. The multi-bosonic representation is taken for 
the first factor P^(x) with a relatively low order n\. This diminishes the storage 
requirements and improves autocorrelations. A better approximation of the fermion 
determinant is achieved by a second polynomial factor P^(x) of order n^. In the 
gauge field update the effect of P^ is taken into account stochastically. Another 
auxiliary polynomial P^\x) is also needed in this stochastic correction. The final 
precision in the approximation is achieved by reweighting the gauge configurations 
when evaluating expectation values. There a fourth polynomial P${x) is used. For 
appropriate algorithms to obtain the necessary optimized polynomial approximations 
see [29~|. A detailed description of the TSMB algorithm can be found in [3D . 

There are two reasons why the TSMB algorithm can overcome the ergodicity 
problem related to the zero of the fermion determinant discussed in section |3.1| . First, 
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the gauge field updates are performed by some large step-size algorithm, in our case 
the Metropolis algorithm. Second, the imperfect approximations near the zero of 
the determinant open a "hole" where the tunnelling between the sectors of differing 
determinant sign is facilitated. The small error in the evaluation of expectation values 
due to the imperfection of the approximation can be removed by the reweighting step 
3~T| . In Fig. ^| we show as a dashed line the evolution of | lndetM| x sign(detM) (for 



the uncorrected M), using a TSMB algorithm and similar simulation parameters to 
the other algorithms in the figure. The algorithm appears to be capable of changing 
the sign of the determinant if anything slightly more effectively than the quenched 
updates. 

A possible procedure to tune the parameters of the TSMB algorithm is as follows: 
since the HMC algorithm is also available, one can determine the smallest (X m in) 
and largest (A max ) eigenvalues of Q 2 on typical gauge configurations from a HMC 
run. The lower and upper limits of the approximation interval [e, A] can be chosen as 
e ~ 0.5A min and A ~ l.h\ max . 

The order of the polynomial for the noisy correction n 2 can be taken to be roughly 
the same as the average number of iterations in the inversions of HMC. The third 
polynomial used in the noisy correction should have an order n 3 which is typically 
10-30% larger than n 2 . (A good test for n 3 is that the noisy correction should ideally 
always accept an unchanged gauge configuration. An acceptance of about 99% is 
sufficient in practice.) After fixing n 2 one can optimise the order of the first polynomial 
by tuning the average acceptance of the noisy correction step. An acceptance of 
about 60-70% turned out to be optimal in most cases. Note that the quality of 
approximation of the fermion determinant is practically independent of ri\ once n 2 is 
fixed. 

An interesting parameter in the optimization of the autocorrelation is the ratio 
of update sweeps performed on the bosonic pseudofermion fields versus gauge fields. 
This depends on the lattice parameters and also on the machine, code optimization, 
compiler etc. In our case we typically found it better to choose two or three times 
as many gauge updates as pseudofermion updates. This differs from the experience 
in super symmetric Yang-Mills theory |3(J where it is better to have relatively more 
boson field updates. 

The final step is to check the approximation of the fermion determinant by calcu- 
lating the reweighting factors on the measured gauge configurations. This can be done 
by determining a few (say, 8 or 16) of the smallest eigenvalues of Q 2 and calculating 
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the reweighting factor explicitly for them, and afterwards multiply by the stochastic 
reweighting factor obtained by a large order (n^) polynomial on the orthogonal sub- 
space. In practice the reweighting does not change the averages if the reweighting 
factors are within a few percent of unity. In runs with large condition numbers A/e, up 
to A/e ~ 10 7 for our parameters, this is not the case and the reweighting is important. 
In fact, in these cases it is not optimal to try to increase so long as the reweighting 
is negligible, since the very high orders required would slow down the updating too 
much. It is better to increase 712 only up to the level that the reweighting factors are 
typically 0(1). The final precision is then achieved by reweighting and this has to be 
done typically only on (more or less) independent configurations. 

Note that if HMC is not available for starting the optimisation one can first 
consider a point with heavy fermions where low order polynomials are sufficient and 
gradually decrease the mass to the point of interest. 

Let us mention a technical point which turns out to be important for dealing with 
large order polynomials, especially in case of a single precision (32 bit) computation. 
The relevant variable for the polynomials is the condition number A/e. The actual 
values of A and e can be reached by rescaling. This enables very large or very small 
numbers appearing in the expansion coefficients of the polynomials to be avoided. It 
turns out that choosing, for instance, A = 4 keeps these numbers within a reasonable 
range. The required rescaling factor can then be included in the recursive evaluation 
of the polynomials. In this way a single precision calculation becomes possible even 
for polynomial orders n = 0(1000). 

4 Results 

4.1 Studies at fi = 

We have performed the bulk of our simulations in this initial study on a 4 3 x 8 
lattice with gauge coupling (3 = 2.0 and quark masses m = 0.1, 0.05 and 0.01. It 
is important first to address the issue raised in subsection |2.4j , namely whether the 
model exhibits confinement and chiral symmetry breaking with these parameters at 
\i = 0, or whether the quarks are already sufficiently light to destroy asymptotic 
freedom and perhaps send the theory into a different phase. 

Fig. |] shows results obtained by HMC simulation at fi = as a function of m. The 
observables monitored, tabulated in Table p], are the chiral condensate (xx) measured 
using a stochastic estimator, the pion mass estimated using a standard cosh fit 
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Figure 4: Results of HMC simulations at (3 — 2.0 on a 4 3 x 8 lattice. 

to all 8 timeslices of the pion propagator, and the average plaquette □ = |Tr (U^), 
which has been rescaled for the convenience of the plot. The data clearly support a 
scenario with lim m _>o(xx( m )) 7^ and m 7r oc yjm in the same limit. This suggests 
that simulations with m > 0.01 fall safely within the regime where chiral symmetry 
is broken; this suffices for our purposes, since chiral symmetry restoration is the main 
physical issue at high density. The increase in the value of the plaquette as m is 
reduced shows nonetheless that color screening due to dynamical fermion effects is 
clearly observable. 

4.2 Autocorrelation Analysis 

Before exploring the phase diagram we performed some long runs in order to de- 
termine the decorrelation time of both algorithms. A crucial question is how the 
simulation effort changes for both HMC and TSMB as one follows the fx axis. To be 
able to compare two different algorithms we need a common unit of measure for the 
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Table 1: 


Results from simulations with jj, 


= 


m 


/ — \ 

\XX) 




□ 


0.10 


1.526(1) 


0.7327(4) 


0.5682(4) 


0.05 


1.453(4) 


0.5298(16) 


0.5805(11) 


0.04 


1.405(6) 


0.4778(14) 


0.5899(15) 


0.03 


1.381(11) 


0.4181(15) 


0.5921(26) 


0.02 


1.307(10) 


0.3482(18) 


0.6004(16) 


0.01 


1.228(8) 


0.2541(25) 


0.6033(14) 



simulation time necessary to obtain two independent configurations. A convenient 
choice is the number of matrix multiplications (appropriately corrected with a factor 
that takes into account that TSMB spends more time in other kinds of operation). 
For HMC the number of matrix multiplications between two successive data takings 
is given by: 

r ptrj 

[(^)(^) + (O x2 x2 (4.1) 

Here T trj is the average trajectory length and dr the elementary step length for the 
hamiltonian deterministic evolution. 1^ is the average number of conjugate gradient 
iterations for a single step in the hamiltonian evolution, and 1^ is the one needed in 
the Metropolis selection step. One factor 2 is there because each conjugate gradient 
iteration implies 2 matrix multiplications, the other because the data are printed out 
every 2 trajectories. 

The corrected number of matrix multiplications to perform a TSMB step is: 

(n 2 + n 3 )xyxF (4.2) 

Here Im is the number of Metropolis iterations in a single TSMB step and F is 
the correction factor mentioned above: in practice F is found to depend on every 
parameter in the simulation, but has typical values 1.2 < F < 2. 

Since we found the plaquette to be the observable with by far the longest auto- 
correlation we concentrated on that for the following analysis. First we considered 
a point at \x = 0. The autocorrelation function for a run with HMC is shown in 
Fig. |5|. We deduce an integrated autocorrelation time of the order of 2 x 10 5 matrix 
multiplications. The corresponding plot for a TSMB run is displayed in Fig. || In this 
case an optimal choice of parameters was found (following the prescription described 
above) to be n\ = 24, n<i = 90, = 110, Im = 12. The integrated autocorrelation 
time turns out to be about 3.2 x 10 5 matrix multiplications. 
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Figure 5: Plaquette autocorrelation function for HMC at /J, — 0.0 

For large values of the chemical potential things get much more difficult for both 
algorithms. For the representative point at m = 0.1, // = 0.4 we cannot give a 
precise determination of the integrated autocorrelations, due to the impracticability 
of accumulating sufficient statistics. However we can see two typical histories in 
Figs. [7| and ^|. The plots look superficially similar; however one can see from the 
number of matrix multiplications that the HMC history is much longer, consisting 
of approximately 15 times more matrix multiplications. This suggests an estimate of 
the slowing down of HMC with respect to TSMB at that point of about one order of 
magnitude. The relatively poor performance of HMC in this region arises from the 
need to reduce dr dramatically (values as small as 0.0002 were needed at the highest 
\i values explored) to maintain reasonable acceptance. These results strongly suggest 
that TSMB may be the algorithm of choice in the high density region, independent 
of the ergodicity considerations discussed in section |[ 

Since the integrated autocorrelation time is not available in the high density re- 
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Figure 6: Plaquette autocorrelation function for TSMB at fi = 0.0 



gion, and also because the determinant sign and reweighting factor need to be taken 
into account, in the following we will quote the jackknife error. For purely gluonic 
observables such as the plaquette this is probably a considerable underestimate of 
the true value, and possibly explains why the mean values we present are not really 
in agreement between the two algorithms. Another factor which may be relevant is 
the lack of ergodicity of HMC, implying that the two algorithms may be exploring 
distinct phase spaces; in subsection [O] we will show that this is the case for some of 
the fermionic observables. 
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Figure 7: Plaquette time history for HMC at /i = 0.4 



4.3 Physics Results from HMC 

Next we report on the results of simulations performed using the HMC algorithm for 
/i ^ 0. We used three distinct quark masses on the 4 3 x 8 lattice at (3 = 2.0, and 
explored values of \x up to and including 0.5 for m = 0.1, \x = 0.45 for m = 0.05, 
and \i = 0.3 for m = 0.01. Our results for the chiral condensate (xx) an d the baryon 
number density n ( 2. 14] ) are shown as a function of \x in Fig. [| 

The results for m = 0.1,0.05 show the condensates remaining unchanged as \x 
increases from zero up to a rather sharply-defined fi , which we identify with the 
onset transition discussed in section |2.1| . At this point the chiral condensate begins 
to fall sharply from its zero-density value, and the baryon density begins to rise 
linearly from zero. The results for m = 0.01, while less clear-cut, are consistent with 
this picture. The lines through the filled points are a straight line fit to the non-zero 
values, the details of which are given in Table ^. The onset value corresponding 
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Figure 8: Plaquette time history for TSMB at fx = 0.4 

to the x-intercept of the fit, coincides quite well with the value of fx for which the 
edge of the eigenvalue spectrum of M(fx) crosses the imaginary axis, as illustrated in 
Fig. § 

The picture is qualitatively very similar to the predictions displayed in figure 4 of 
||. We can make the comparison more quantitative by replotting the data in terms 
of rescaled variables x = 2fi/m 7T0 , y = (xx) / (xx)o, an d n = r nnon/8m(xx)oi where 
the subscript denotes values at zero chemical potential. The predictions from %PT 
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Figure 9: (xx) (open symbols) and n (filled symbols) vs. /i, for masses m = 0.1 
(circles), m = 0.05 (squares) and m = 0.01 (diamonds). 



are then that all data should fall on the lines 




n 



x (l _l 



X < 1, 



X > 1. 



(4.3) 



Using the data of Fig. ^ and Table [l] we plot y vs. x in Fig. 10 and n vs. x in 
Fig. |TT|. The data collapse very nicely onto a universal curve, corresponding quite 
closely to the prediction flOp . The systematic departures from the theoretical curves, 
downwards for the condensate data and upwards for the baryon density, may well be 
explicable by higher order corrections in x?T. In Table [2] we also list predictions 
for the slope and x-intercept of the linear fit to the baryon density, using a linearised 
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Table 2: Linear fits to in the vicinity of the onset transition 



m 


x-Intercept 


x-Intercept 


Slope 


Slope 




(XPT) 


(measured) 


(XPT) 


(measured) 


0.10 


0.3664(2) 


0.356(8) 


4.548(6) 


3.96(6) 


0.05 


0.2649(8) 


0.249(5) 


4.141(28) 


2.94(4) 


0.01 


0.1271(13) 


0.117(11) 


3.044(69) 


1.84(8) 



approximation to (|4.3|) : 



n - 16^%^ L _ ^ ) , (4 . 4) 



m 



7r0 



The match between theory and measurement is excellent for the intercept and strongly 
supports the identification \i a ~ m n /2; that for the slope less so. In any case, the 
agreement between theory and measurement is remarkable for data taken on such a 
small lattice, relatively far from the continuum limit, with quark masses ranging over 
an order of magnitude. 

It seems reasonable to deduce that the high density phase for /i > /i is superfluid, 



characterised by a non-zero diquark condensate of the form qqs ( 2.16|) , similar to that 
observed in lattice simulations of Two Color QCD with fundamental fermions [[15 . 
Work to establish this by direct measurement is in progress. Whilst the quantitative 
agreement between our results and the theoretical predictions of |9|] is gratifying, 
it also contradicts the symmetry-based arguments of section |2.1| that there are no 
baryonic Goldstones for N — 1 staggered flavor, and no gauge-invariant local diquark 
condensate. We believe that this is because the HMC simulations fail to take into 
account of the determinant sign (or indeed even to change it) ie. that simulations with 
functional weight |detM| yield broadly similar results to those with weig ht det 2 M; 
the premature onset at /i G = is therefore a direct manifestation of the sign 

and/or ergodicity problems. 

Next, we investigate function of //. Recall that for // ^ 0, the pion 

timeslice propagator is defined as 



necessitating two inversions of M. We implemented fl4.5|) using a source site and color 
chosen at random and summing over sink color, performing a simple cosh fit to the 
resulting G n over all 8 timeslices. The fits were stable in the low-density phase, but 
above the onset transition G n became markedly noisier and the fit less convincing, 
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Figure 10: Chiral condensate vs. chemical potential using the rescaled variables of 
eq. (|4.3p . The symbols are the same as those of Fig. |^ 

suggesting that perhaps a different functional form is more suitable. Our results for 
m n are shown in Fig. [E| For fi < fi is constant to quite high precision; for 
\l > // it begins to rise. Qualitatively similar behaviour has been observed in the 
Gross-Neveu model WA \ however in the present case we also have the theoretical 
treatment of ||, which predicts that in the high density phase the state with the 
quantum numbers of the pion has mass m v = 2fi. Whilst the large errors in the 
dense phase preclude a precise comparison, the clustering of the points just above the 
line m-x = 2fi is striking. 

The issue of whether chiral symmetry is restored in the dense phase, ie. whether 
lim m ^ (xx( m )) 7^ 0, is complicated by the sensitivity of /i G to m oc m^. Eq. (|4~3| ) 
suggests that (xx) should decrease as fi~ 2 for \i > approaching zero only asymp- 
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Figure 11: Baryon density vs. chemical potential using the rescaled variables of eq. 

Q. 

totically as \i —>■ oo. This follows from the idea that the chiral condensate is gradually 
rotated into a diquark condensate as /i increases. Fig. [It] shows the first hints of this 
behaviour. We can, however, take a more pragmatic (as well as more physical) ap- 
proach and plot (xx{ m )) a t fixed n rather than at fixed \i. This has been done for 
n = 0.1 0.2 and 0.3 using the linear fit to n(/i) of Table |2] and a simple-minded linear 
interpolation of the chiral condensate data. The result is shown in Fig. |13[ Clearly 
data from more values of m would be needed to make a definitive statement, but 
there is a suggestion that chiral symmetry is not completely restored, particularly at 
the lowest density n = 0.1. 

Finally we turn to the effect of the chemical potential on the gauge fields. Since 
this can only be communicated via fermion loops, any effect we see can be ascribed 
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Figure 12: vs. //, for masses m = 0.1 (circles), m = 0.05 (squares) and m = 0.01 
(diamonds). Also shown is the line m n = 2/i. 



with certainty to dynamical fermions (we cannot exclude the possibility that all our 
other observations could have been made for a fraction of the cost in the quenched 
approximation). Gluonic observables, however, are also much more prone to auto- 
correlations as described in section |47|, particularly as the quark mass is reduced. 
Systematic changes with /i are therefore quite difficult to observe. In this initial 
HMC study we have only measured the average plaquette; the results are shown in 



Fig. [14]. The data for m = 0.1, 0.05 show the plaquette remaining roughly constant 
for jj, < fi D , before beginning to decrease. The m = 0.01 data are consistent with this 
picture within admittedly large errors. We interpret it as follows: for temperature 
T = 0, all values of /i < fi are physically equivalent corresponding to the same 
physical state, namely the vacuum. We only expect an effect on gluonic observables 
in the presence of matter, ie. for \x > \x Q . To the extent that the results are constant 
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Figure 13: (xx) vs - m f° r various fixed baryon densities. 

for /i < fi Q we can be confident that our simulation has an effective T ~ 0. The 
decrease in the plaquette for ji > \i may be due to the decrease in the number 
of virtual quark-anti-quark pairs which may form due to the Exclusion Principle 
- an effect known as Pauli Blocking. This results in a decrease of screening via 
vacuum polarisation, and hence an effective renormalisation of the gauge coupling 
(5 and consequent decrease of the plaquette. In the large-/! limit the lattice should 
become saturated with one quark of each color per site, and the plaquette assume its 
quenched value; this has been verified at strong gauge coupling . 
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Figure 14: Average plaquette □ vs. fi for m 
(diamonds). 



0.1 (circles), 0.05 (squares) and 0.01 



4.4 Physics Results from TSMB 

Each TSMB simulation is characterised by a vector rij specifying the polynomial 
orders at each stage, as described in section j3l2[ . For each configuration generated, 
a reweighting factor r and the sign of detM must be determined. On the relatively 
small lattices considered here, it is possible to compute detM directly using standard 
numerical methods; on larger lattices one can use the 'spectral flow method' [3(| as a 
function of m and/or \x. The expectation value of an observable O is then determined 
by the ratio 

(Oxrx sign) 



(O) 



(r x sign) 



(4.6) 
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Here we present results from runs on a 4 3 x 8 lattice with (3 = 2.0, m = 0.1 at 
three values of chemical potential: /x = 0.0 with polynomials n« = (64,250,300); 
fi = 0.36 with both m = (64, 300, 350, 500) and n, = (80, 300, 350, 500); and /i = 0.4 
with rii = (160, 1000, 1100, 1200) (reweighting was not performed at /i = 0.0). Note 
that the polynomial orders required increase with /i. The point at \i = was chosen 
to enable the TSMB algorithm to be tested against HMC, since both should yield 
identical results. The /i ^ points were chosen so as to have one value just past 
the HMC onset transition, where the edge of the eigenvalue distribution just overlaps 
the line ReA = and a small percentage of negative determinant configurations are 
expected, so that hopefully the sign problem is not too severe, and one value fairly 
deep in the high density phase. 

Our results for the standard observables, together with the corresponding HMC 
results, are summarised in Table [4.4| . For TSMB at /i ^ we also include observables 
determined separately in each sign sector, defined by (0)± = (O x r)±/(r)±. The 
observables at \x = 0.36 and 0.4 result from runs on 32 separate configurations with a 
few thousand update cycles each. At \i = 0.0 the results quoted for (xx) an d n come 
from a run on 128 configurations with ~1500 update cycles. Autocorrelation studies 
reveal that these may not be fully thermalised or decorrelated, which may account 
for the slight discrepancies between TSMB and HMC. The plaquette at ji = 0.0, by 
far the slowest observable to decorrelate, is taken from a long run of ~25000 update 
cycles on a single configuration with the same Wj, with measured autocorrelations 
taken into account, but ignoring reweighting (reweighting, including sign(detM) has 
not been observed to have a significant effect on the plaquette). For \i = 0.0 and 
0.4 the second number quoted results from long runs with 40000 update cycles at 
Hi = (24, 90, 120) and 19000 at m = (140, 1000, 1100) — these runs have been used to 
obtain Figs. |6| and || respectively. Even so, at \i — 0.4 the long autocorrelation time 
implies that the errors are likely to be underestimated, which perhaps explains the 
large discrepancy between the plaquette value including the measured autocorrelation 
and the other plaquette values (note also that the HMC plaquette value at /i = 0.36 
appears to lie outside the trend of Fig. |i~4l) . 

At n = 0.0 the plaquette results in the second row provide reassurance that the 
two algorithms agree, and that we have some degree of control over errors. Only once 
/i > Ho (as determined by HMC) do we anticipate the two should differ as a result of 



HMC's failure to explore the negative determinant sector. In Figs. [15] and we plot 
distributions of r x sign obtained at \i = 0.36, 0.4 respectively. For fi = 0.36 there is 
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fi 




TSMB 




HMC 






(0) 


(0) + 


(0)- 




(xx) 


0.0 


1.510(4) 






1.526(1) 




0.36 


1.562(18) 


1.538(15) 


1.240(58) 


1.485(9) 




0.4 


1.26(15) 


1.24(4) 


1.23(4) 


1.253(10) 


n 


0.0 


-0.0018(26) 






-0.0002(3) 




0.36 


-0.005(14) 


0.015(12) 


0.263(52) 


0.0172(28) 




0.4 


0.09(13) 


0.17(3) 


0.21(3) 


0.1667(90) 


□ 


0.0 


0.5765(58) 






0.5682(4) 






0.5667(17) 










0.36 


0.5541(12) 


0.5542(11) 


0.5551(18) 


0.5729(40) 




0.4 


0.588(7) 


0.5791(14) 


0.5746(15) 


0.5612(30) 






0.5523(46) 









Table 3: A comparison of results between TSMB and HMC 



a small fraction of negative determinant configurations; for fi = 0.4 the distribution 
is much more symmetrical between + and -, although the precise shape is found 
to be sensitive to the choice of rii, and ri2 in particular. The severity of the sign 
problem in numerical simulations is usually expressed in terms of the average sign 
(sign); generically this decreases towards zero as the system volume increases until 
its relative error becomes so large that accurate estimates of (O) are impracticable. 
For TSMB the corresponding quantity is not uniquely defined; eg. for the data of 
Fig. [15| we could specify the fraction of negative determinant configurations (12%), 
(r x sign)/(r) = 0.910, or ((r) + - (r)_)/«r) + + (r)_) = 0.491. 

A signal for physical effects associated with the inclusion of the sign of the deter- 
minant in the functional measure following ( f4.6| ) is that (0)+ ^ The centre 
columns of table [4.4| reveal evidence for such an effect in the fermionic obsevables 
for fi = 0.36. As a result, the full average over both sectors for (xx) is significantly 
greater than the HMC result, and perhaps even consistent with (xx)o- More spec- 
tacularly, the average baryon density n is consistent with zero. These observations 
imply that at this value of fi the system is still in the low density phase, and hence 
Ho true > Ho hmc- This is in accord with our symmetry-based arguments that for 
N — 1 flavors of adjoint staggered fermion there are no baryonic Goldstones and 
hence no early onset. Note that similar effects seen at fi — 0.4 are statistically far 
less significant, and should at this stage be considered preliminary. 

We have seen that TSMB simulations give evidence for a delayed onset, giving us 
confidence that the algorithm correctly samples the two sign sectors and thus correctly 
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Figure 15: Histogram of r x sign for = 0.36 with n; = (64, 300, 350, 500). 
describes a single flavor. This is the principal result of our initial TSMB studies. 

5 Conclusions 

To our knowledge this has been the first study of Two Color QCD with adjoint quarks, 
and the first TSMB study both to use staggered fermions, and to set /i ^ 0. The 
highlights of this work are the following: 



• We have outlined the global symmetries and anticpated patterns of symmetry 
breaking for SU(2) lattice gauge theory with staggered fermions and a non-zero 
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Figure 16: Histogram of r x sign for /i = 0.40 with = (200, 1000, 1100, 1200). 



chemical potential in both fundamental and adjoint representations. The case 
of N = 1 adjoint flavor seems especially interesting, being the most 'QCD-like'. 

• We have studied the model using both conventional hybrid Monte Carlo (HMC) 
and Two-Step Multi-Bosonic (TSMB) algorithms. The HMC simulations slow 
down dramatically, both in terms of number of matrix inversions required, and 
in terms of autocorrelation, once the eigenvalue distribution includes some with 
negative real part, which begins to occur once baryon density n > 0. We 
have confirmed the existence of isolated real negative eigenvalues at large /i and 
hence a sign problem. The HMC algorithm is unable to change the sign of the 
determinant, and hence is not ergodic in this region. The TSMB algorithm, 
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by reason of the approximate way in which it treats small eigenvalues, is able 
to change the determinant sign — our first studies also indicate it decorrelates 
more effectively than HMC in the dense phase. 

Simulations using HMC over a range of m and /x show good quantitative agree- 
ment with a chiral perturbation theory treatment — data spanning an order 
of magnitude in quark mass collapsing onto a simple universal curve whose pa- 
rameters are completely determined by physical quantities measured at /x = 0. 
In particular a second order phase transition (see Fig. [TTJ) at /x — m n /2 sepa- 



rates the vacuum from a phase with n > 0, which can be interpreted as a fluid 
of mutually repelling diquark bosons 0. The xPT analysis should only hold 
for the global symmetries of the model with N > 2 flavors, implying that the 
failure to explore the negative determinant sector of the model results in the 
'wrong' physics. 

• Measurements made using TSMB simulations which take the determinant sign 
into account give evidence for a delayed onset, ie. n remains consistent with 
zero even for /x > m n /2. This is in accord with the symmetry breaking pattern 
anticipated for N = 1. 

In the future we plan to extend our simulations of the dense phase using both 
algorithms, with the following goals: 

• We wish to examine the spectrum of the model at /x = in greater detail, 
in order to establish the masses of baryonic states such as vector diquarks, 
qg fermions, etc, which may be associated with further thresholds as they are 
induced into the ground state as /x is raised. We expect the %PT treatment, 
which does not include such states, to cease to be accurate at some point. This 
may be signalled either by the breakdown of universality of data with different 
m, or by further phase transitions. 

• We wish to probe larger values of /x, up to the saturation point n ~ 3 quarks 
per lattice site, to see whether any new phases emerge. Perhaps at some point 
a sharp Fermi surface will appear. 

• We wish to explore gluodynamics in the dense medium []10| |, by studying (in 



order of sophistication) Pauli Blocking, the static quark potential, the gluon 
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propagator, and the spatial and size distribution of instantons. These studies 
will require a sizeable increase in lattice volume. 

• We wish to explore signals for diquark condensation in the dense phase, both in 
superfluid and superconducting channels. This latter exotic phenomenon will 
require good quantitiative control over the determinant sign fluctuations, and 
hence high statistics (though probably not a large volume). 
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